module powermethod

contains

subroutine powermethodG(m,n, mt, phi, psi, funcFS, funcQ, lsolver, ccpsi, k, err, dr, iter)
    use LU_tdiag
    use ifport
    implicit none
    include "xsec.FI"
    include "files.FI"
    interface 
        subroutine funcFS
        end subroutine
        subroutine funcQ(g, k, q)
            integer, intent(in) :: g
            real(8), intent(in) :: k
            real(8), intent(out), dimension(:) :: q
        end subroutine
        subroutine lsolver(sm, mt, rhs, sol)
            integer, intent(in) :: sm
            real(8),dimension(3,sm), intent(in) :: mt
            real(8),dimension(sm), intent(in) :: rhs
            real(8),dimension(sm), intent(out) :: sol
        end subroutine
        
    end interface

    integer, intent(in) :: m,n
    real(8), dimension(3,m,ng), intent(in):: mt
    real(8), intent(in):: ccpsi
    
    real(8), dimension(m,ng), intent(inout):: phi
    real(8), dimension(m), intent(inout):: psi
    real(8), dimension(m):: psi_prv
    real(8), dimension(m):: q
    real(8), intent(inout):: k
    real(8), intent(out) ::  err, dr
    integer, intent(out) :: iter
    
    integer i, j, g, gp
    real(8) torErr, upper, lower
    real(8) :: err_prv
    
    iter   = 0
    err    = 1.
    dr = 1.
!    do g=1, ng
!        do i=1, m
!            phi(i,g) = drand(0)
!        end do
!    end do
!    phi=1.
    
    ! make fission source
    call funcFS
    
    do while(err.gt.ccpsi)
        err_prv = err
        iter = iter + 1
        do g=1, ng
            ! make total source
            call funcQ(g,k,q)
            ! solve 
            call lsolver (m, mt(:,:,g), q, phi(:,g))
        end do
            
        ! make fission source from updated flux
        psi_prv = psi
        call funcFS
                
        ! evaluate k                                            
        upper = 0.
        lower = 0.
        do i=1, m                                                
            upper = upper + psi(i)*psi(i)
            lower = lower + psi(i)*psi_prv(i)
        end do
        k = k*upper/lower
        
        ! evaluate err 
        err = 0.
        lower = sqrt(upper)
        upper = 0
        do i=1, m
            upper = upper + (psi(i)- psi_prv(i))*(psi(i)- psi_prv(i))
        end do
        upper = sqrt(upper)
        err = upper/lower
        if (iter==1000) exit

    end do
    dr = err/err_prv
end subroutine

subroutine powermethodF(m,n, mt, phi, psi, funcFS, funcQ, lsolver,ccpsi, k, err, dr, iter)
    use LU_tdiag
    use ifport
    implicit none
    include "xsec.FI"
    include "files.FI"
    interface 
        subroutine funcFS
        end subroutine
        subroutine funcQ(k, q)
            real(8), intent(in) :: k
            real(8), intent(out), dimension(:) :: q
        end subroutine
        subroutine lsolver(size_, mtx_, rhs_, sol_)
            integer, intent(in) :: size_
            real(8), intent(in), dimension(:,:) :: mtx_
            real(8), intent(in), dimension(:) :: rhs_, sol_
        end subroutine
        
    end interface

    integer, intent(in) :: m,n
    real(8), dimension(n,m), intent(in):: mt
    real(8), intent(in):: ccpsi
    
    real(8), dimension(m,ng), intent(inout):: phi
    real(8), dimension(m), intent(inout):: psi
    real(8), dimension(m):: psi_prv
    real(8), dimension(m*ng) :: q, phiar
    real(8), intent(inout):: k
    real(8), intent(out) ::  err, dr
    integer, intent(out) :: iter
    
    integer i, j, g, gp
    real(8) torErr, upper, lower
    real(8) :: err_prv
    
    iter   = 0
    err    = 1.
    dr = 1.
!    do g=1, ng
!        do i=1, m
!            phi(i,g) = drand(0)
!        end do
!    end do
!    phi=1.
    
    ! make fission source
    call funcFS
    
    do while (err.gt.ccpsi)
        err_prv = err
        iter = iter + 1
        ! make total source
        call funcQ(k,q)
        ! solve 
        call lsolver(m*ng, mt, q, phiar)
        do g=1, ng
            phi(:,g) = phiar((g-1)*m+1 : g*m)
        end do
        ! make fission source from updated flux
        psi_prv = psi
        call funcFS
                
        ! evaluate k                                            
        upper = 0.
        lower = 0.
        do i=1, m                                                
            upper = upper + psi(i)*psi(i)
            lower = lower + psi(i)*psi_prv(i)
        end do
        k = k*upper/lower
        
        ! evaluate err 
        err = 0.
        lower = sqrt(upper)
        upper = 0
        do i=1, m
            upper = upper + (psi(i)- psi_prv(i))*(psi(i)- psi_prv(i))
        end do
        upper = sqrt(upper)
        err = upper/lower
        if (iter==30000) exit

    end do
    dr = err/err_prv
end subroutine    

subroutine powermethodW(m, n, mt, phi, psi, funcFS, funcQ, shift, lsolver, ccpsi, k, err, dr, iter)
    use ifport
    implicit none
    include "xsec.FI"
    include "files.FI"
    interface 
        subroutine funcFS
        end subroutine
        subroutine funcQ(k,q)
            real(8), intent(in) :: k
            real(8), dimension(:), intent(out) :: q
        end subroutine
        subroutine shift(ke)
            real(8), intent(in) :: ke
        end subroutine
        subroutine lsolver(sm, mt, rhs, sol)
            integer, intent(in) :: sm
            real(8),dimension(sm,sm), intent(in) :: mt
            real(8),dimension(sm), intent(in) :: rhs
            real(8),dimension(sm), intent(out) :: sol
        end subroutine
    end interface

    integer, intent(in) :: m,n
    real(8), dimension(n, m*ng), intent(in):: mt
    real(8), intent(in):: ccpsi
    real(8), dimension(m,ng), intent(inout):: phi
    real(8), dimension(m), intent(inout):: psi
    real(8), dimension(m):: psi_prv
    real(8), dimension(m*ng):: q, phiar
    real(8), intent(inout):: k
    real(8), intent(out) ::  err, dr
    integer, intent(out) :: iter
    
    integer i, j, g, gp
    real(8) torErr, upper, lower
    real(8) :: err_prv, ks, ke, kd
    
    iter   = 0
    err    = 1.
    dr = 1.
    phi=1.
    kd =0.01
    
    ! make fission source
    call funcFS
    phiar = 0;
    do while (err.gt.ccpsi)
        err_prv = err
        iter = iter + 1
        ! make total source
        ke = k + shiftk(iter,kd)
        ks = 1._8/(1._8/k-1._8/ke)
        call shift(ke)
        call funcQ(ks,q)
        ! solve 
        call lsolver(m*ng, mt(:,:), q, phiar)
        do g=1, ng
            phi(:,g) = phiar((g-1)*m+1 : g*m)
        end do
            
        ! make fission source from updated flux
        psi_prv = psi
        call funcFS
                
        ! evaluate k                                            
        upper = 0.
        lower = 0.
        do i=1, m                                                
            upper = upper + psi(i)*psi(i)
            lower = lower + psi(i)*psi_prv(i)
        end do
        ks = ks*upper/lower
        k = (ks*ke)/(ks+ke)
                
        ! evaluate err 
        err = 0.
        lower = sqrt(upper)
        upper = 0
        do i=1, m
            upper = upper + (psi(i)- psi_prv(i))*(psi(i)- psi_prv(i))
        end do
        upper = sqrt(upper)
        err = upper/lower
        if (iter==100) exit

    end do
    dr = err/err_prv  
end subroutine

subroutine powermethodWT(m, n, mt, phi, psi, funcFS, funcQ, shift, lsolver, ccpsi, k, err, dr, iter)
    use ifport
    implicit none
    include "xsec.FI"
    include "files.FI"
    interface 
        subroutine funcFS
        end subroutine
        subroutine funcQ(k,q)
            real(8), intent(in) :: k
            real(8), dimension(:), intent(out) :: q
        end subroutine
        subroutine shift(ke)
            real(8), intent(in) :: ke
        end subroutine
        subroutine lsolver(sm, mt, rhs, sol)
            integer, intent(in) :: sm
            real(8),dimension(3,sm), intent(in) :: mt
            real(8),dimension(sm), intent(in) :: rhs
            real(8),dimension(sm), intent(out) :: sol
        end subroutine
    end interface

    integer, intent(in) :: m,n
    real(8), dimension(n, m*ng), intent(in):: mt
    real(8), intent(in):: ccpsi
    real(8), dimension(m,ng), intent(inout):: phi
    real(8), dimension(m), intent(inout):: psi
    real(8), dimension(m):: psi_prv
    real(8), dimension(m*ng):: q, phiar
    real(8), intent(inout):: k
    real(8), intent(out) ::  err, dr
    integer, intent(out) :: iter
    
    integer i, j, g, gp
    real(8) torErr, upper, lower
    real(8) :: err_prv, ks, ke, kd
    
    iter   = 0
    err    = 1.
    dr = 1.
    phi=1.
    kd =0.01
    
    ! make fission source
    call funcFS
    phiar = 0;
    do while (err.gt.ccpsi)
        err_prv = err
        iter = iter + 1
        ! make total source
        ke = k + shiftk(iter,kd)
        ks = 1._8/(1._8/k-1._8/ke)
        call shift(ke)
        call funcQ(ks,q)
        ! solve 
        call lsolver(m*ng, mt(:,:), q, phiar)
        do g=1, ng
            phi(:,g) = phiar((g-1)*m+1 : g*m)
        end do
            
        ! make fission source from updated flux
        psi_prv = psi
        call funcFS
                
        ! evaluate k                                            
        upper = 0.
        lower = 0.
        do i=1, m                                                
            upper = upper + psi(i)*psi(i)
            lower = lower + psi(i)*psi_prv(i)
        end do
        ks = ks*upper/lower
        k = (ks*ke)/(ks+ke)
                
        ! evaluate err 
        err = 0.
        lower = sqrt(upper)
        upper = 0
        do i=1, m
            upper = upper + (psi(i)- psi_prv(i))*(psi(i)- psi_prv(i))
        end do
        upper = sqrt(upper)
        err = upper/lower
        if (iter==1000) exit

    end do
    dr = err/err_prv  
end subroutine    

subroutine powermethod_2nd (size, mt, funcFS, funcQ, ccpsi, phi, psi, k, err, dr, iter)
    use LU_tdiag
    use ifport
    implicit none
    include "xsec.FI"
    include "files.FI"
    interface 
        subroutine funcFS
        end subroutine
        subroutine funcQ(g,k, q)
            integer, intent(in) :: g
            real(8), intent(in) :: k
            real(8),intent(out),dimension(:) ::q
        end subroutine
    end interface

!    external funcFS ! fission source
!    external funcQ  ! total source 

    integer, intent(in) :: size
    real(8), dimension(3,size,ng), intent(in):: mt
    real(8), intent(in):: ccpsi
    
    real(8), dimension(size,ng), intent(inout):: phi
    real(8), dimension(size), intent(inout):: psi
    real(8), dimension(size):: q
    real(8), intent(inout):: k
    real(8), intent(out) ::  err, dr
    integer, intent(out) :: iter
    
    real(8), dimension(size):: psi_prv
    integer i, j, g, gp
    real(8) torErr, upper, lower
    real(8) :: err_prv
    
    real(8) :: phi1st(size, ng)
    real(8) :: ci, l2nm1, avg
    
    phi1st = phi
    do g=1, ng
        do i=1, size/2
            avg = (phi1st(i,g) + phi1st(size+1-i, g))/2
            phi1st(i,g) = avg
            phi1st(size+1-i, g) = avg
        end do
    end do
!    phi = 1.
    do g=1, ng
        do i=1, size
            phi(i,g) = drand(0)
        end do
    end do
    l2nm1 = inner_prod(ng, size, phi1st, phi1st)
    ci = inner_prod(ng, size, phi, phi1st)/l2nm1
    phi = phi - ci * phi1st
    call even
    
    iter   = 0
    err    = 1.
    dr = 1.
    
    ! make fission source
    call funcFS
    
    do while (err.gt.ccpsi)
        err_prv = err
        iter = iter + 1
        do g=1, ng
            ! make total source
            call funcQ(g,k,q)
            ! solve 
            call solve_LU (size, mt(:,:,g), q, phi(:,g))
        end do
        ci = inner_prod(ng, size, phi, phi1st)/l2nm1
        phi = phi - ci * phi1st
        call even            
        
        ! make fission source from updated flux
        psi_prv = psi
        call funcFS
                
        ! evaluate k                                            
        upper = 0.
        lower = 0.
        do i=1, size                                                
            upper = upper + psi(i)*psi(i)
            lower = lower + psi(i)*psi_prv(i)
        end do
        k = k*upper/lower
        
        ! evaluate err 
        err = 0.
        lower = sqrt(upper)
        upper = 0
        do i=1, size
            upper = upper + (psi(i)- psi_prv(i))*(psi(i)- psi_prv(i))
        end do
        upper = sqrt(upper)
        err = upper/lower
        if (iter==30000) exit

    end do
    dr = err/err_prv  
    
    
contains

function inner_prod (dim1, dim2, v1, v2) result(ret)
    integer, intent(in) :: dim1, dim2
    real(8), dimension(dim2, dim1) :: v1,v2
    integer:: i, j
    real(8) :: ret
    ret=0._8
    do i=1, dim1
        do j=1, dim2
            ret = ret + v1(j,i)*v2(j,i) 
        enddo
    enddo
    
end function

subroutine even
!    do g=1, ng
!        do i=1, size/2
!            avg = (phi(i,g) + phi(size+1-i, g))/2
!            phi(i,g) = avg
!            phi(size+1-i, g) = avg
!        end do
!    end do
end subroutine
        
end subroutine
        
        
!subroutine powermethod_wldt(size, mt, funcFS, funcQ, matrix_shift, ccpsi, phi, psi, k, err, dr, iter)
!    use LU_general
!    use ifport
!    implicit none
!    include "xsec.FI"
!    include "files.FI"
!    interface 
!        subroutine funcFS
!        end subroutine
!        subroutine funcQ(k,q)
!            real(8), intent(in) :: k
!            real(8), dimension(:), intent(in) :: q
!        end subroutine
!        subroutine matrix_shift(ke)
!            real(8), intent(in) :: ke
!        end subroutine
!    end interface
!
!    integer, intent(in) :: size
!    real(8), dimension(size*ng, size*ng), intent(in):: mt
!    real(8), intent(in):: ccpsi
!    
!    real(8), dimension(size,ng), intent(inout):: phi
!    real(8), dimension(size), intent(inout):: psi
!    real(8), dimension(size*ng):: q, phiar
!    real(8), intent(inout):: k
!    real(8), intent(out) ::  err, dr
!    integer, intent(out) :: iter
!    
!    real(8), dimension(size):: psi_prv
!    integer i, j, g, gp
!    real(8) torErr, upper, lower
!    real(8) :: err_prv, ks, ke, kd
!    
!    iter   = 0
!    err    = 1.
!    dr = 1.
!!    phi=1.
!    kd =0.01
!    
!    ! make fission source
!    call funcFS
!    phiar = 0;
!    do while (err.gt.ccpsi)
!        err_prv = err
!        iter = iter + 1
!        ! make total source
!        ke = k + kd
!        ks = 1._8/(1._8/k-1._8/ke)
!        call matrix_shift(ke)
!        call funcQ(ks,q)
!        ! solve 
!        call solve_LU (size*ng, mt(:,:), q, phiar)
!        do g=1, ng
!            phi(:,g) = phiar((g-1)*size+1 : g*size)
!        end do
!            
!        ! make fission source from updated flux
!        psi_prv = psi
!        call funcFS
!                
!        ! evaluate k                                            
!        upper = 0.
!        lower = 0.
!        do i=1, size                                                
!            upper = upper + psi(i)*psi(i)
!            lower = lower + psi(i)*psi_prv(i)
!        end do
!        ks = ks*upper/lower
!        k = (ks*ke)/(ks+ke)
!                
!        ! evaluate err 
!        err = 0.
!        lower = sqrt(upper)
!        upper = 0
!        do i=1, size
!            upper = upper + (psi(i)- psi_prv(i))*(psi(i)- psi_prv(i))
!        end do
!        upper = sqrt(upper)
!        err = upper/lower
!        if (iter==30000) exit
!
!    end do
!    dr = err/err_prv  
!end subroutine
!
!
!subroutine powermethod_wldt2(size, mt, funcFS, funcQ, shift, solver, ccpsi, phi, psi, k, err, dr, iter)
!    use ifport
!    implicit none
!    include "xsec.FI"
!    include "files.FI"
!    interface 
!        subroutine funcFS
!        end subroutine
!        subroutine funcQ(k,q)
!            real(8), intent(in) :: k
!            real(8), dimension(:), intent(in) :: q
!        end subroutine
!        subroutine shift(ke)
!            real(8), intent(in) :: ke
!        end subroutine
!        subroutine solver(sm, mt, rhs, sol)
!            integer, intent(in) :: sm
!            real(8),dimension(3,sm), intent(in) :: mt
!            real(8),dimension(sm), intent(in) :: rhs
!            real(8),dimension(sm), intent(out) :: sol
!        end subroutine
!    end interface
!
!    integer, intent(in) :: size
!    real(8), dimension(3, size*ng), intent(in):: mt
!    real(8), intent(in):: ccpsi
!    
!    real(8), dimension(size,ng), intent(inout):: phi
!    real(8), dimension(size), intent(inout):: psi
!    real(8), dimension(size*ng):: q, phiar
!    real(8), intent(inout):: k
!    real(8), intent(out) ::  err, dr
!    integer, intent(out) :: iter
!    
!    real(8), dimension(size):: psi_prv
!    integer i, j, g, gp
!    real(8) torErr, upper, lower
!    real(8) :: err_prv, ks, ke, kd
!    
!    iter   = 0
!    err    = 1.
!    dr = 1.
!!    phi=1.
!    kd =0.01
!    
!    ! make fission source
!    call funcFS
!    phiar = 0;
!    do while (err.gt.ccpsi)
!        err_prv = err
!        iter = iter + 1
!        ! make total source
!        ke = k + kd
!        ks = 1._8/(1._8/k-1._8/ke)
!        call shift(ke)
!        call funcQ(ks,q)
!        ! solve 
!        call solver (size*ng, mt(:,:), q, phiar)
!        do g=1, ng
!            phi(:,g) = phiar((g-1)*size+1 : g*size)
!        end do
!            
!        ! make fission source from updated flux
!        psi_prv = psi
!        call funcFS
!                
!        ! evaluate k                                            
!        upper = 0.
!        lower = 0.
!        do i=1, size                                                
!            upper = upper + psi(i)*psi(i)
!            lower = lower + psi(i)*psi_prv(i)
!        end do
!        ks = ks*upper/lower
!        k = (ks*ke)/(ks+ke)
!                
!        ! evaluate err 
!        err = 0.
!        lower = sqrt(upper)
!        upper = 0
!        do i=1, size
!            upper = upper + (psi(i)- psi_prv(i))*(psi(i)- psi_prv(i))
!        end do
!        upper = sqrt(upper)
!        err = upper/lower
!        if (iter==30000) exit
!
!    end do
!    dr = err/err_prv  
!end subroutine
function shiftk(iter, delk) result (ret)
    real(8) :: SHIFT_GUIDe(17)
    data SHIFT_GUIDE /1.0_8, 0.75_8, 0.50_8, 0.35_8, 0.25_8, &
                    0.20_8, 0.15_8, 0.125_8, 0.100_8, 0.075_8, 0.050_8, &
                    0.035_8, 0.025_8, 0.020_8, 0.015_8, 0.0125_8, 0.0100_8/
!    data SHIFT_GUIDE /0.5_8, 0.25_8, 0.1_8, 0.05_8, 0.01_8/
    integer, intent(in) :: iter
    real(8), intent(in) :: delk
    real(8) :: ret
!    ret = 1000;
!    return
    ret = delk
    if (iter .le. 17) then
        if (ret<SHIFT_GUIDE(iter)) ret = SHIFT_GUIDE(iter) 
    endif
!    ret = ret * 100
end function

end module